A global flow scheme for MD in \brahms is given in Figure~\ref{fig:mdGlobal}.
\begin{figure}
\begin{center}
\addtolength{\fboxsep}{0.5cm}
\begin{shadowenv}[12cm]
{\large \bf THE BRAHMS MD ALGORITHM}
\rule{\textwidth}{2pt} \\
{\bf 1. Input initial conditions}\\[2ex]
Force-field and simulation parameters\\
Positions and velocities of all sites in the system \\
$\Downarrow$\\
\rule{\textwidth}{1pt}\\
{\bf repeat 2,3,4} for the required number of steps:\\
\rule{\textwidth}{1pt}\\
{\bf 2. Compute forces and torques} \\[1ex]
Energies and pressure are also computed and averaged at a certain frequency. \\   
$\Downarrow$\\
{\bf 3. Update configuration - Integration} \\[1ex]
The movement of the particles is simulated by numerically solving Newton's and Euler's equations of motion \\[1ex]
{\bf 4.} {\bf Output step} \\
write positions, velocities, energies, temperature, pressure, etc. \\
\end{shadowenv}
\caption[The molecular dynamics algorithm]{The global molecular dynamics algorithm as implemented in \brahms.}
\label{fig:mdGlobal}
\end{center}
\end{figure}

 The integration of the equations of motion can be done by using the simplest (possibly) of numerical techniques: the Verlet scheme~\cite{verlet67}. In general, Verlet-like integrators  are time-reversible and symplectic\footnote{Symplecticity implies that volume in phase space is conserved~\cite{allen04}.}: time reversibility and symplecticity have been shown to produce a favourable propagation of errors, hence better numerical methods, especially when interested in long-term behaviour~\cite{mclac05a}.
All the algorithms presented in the following sections are explicit, i.e. not requiring any iteration. 
It should be pointed out that the speed of the integrator itself is not very relevant, because the fraction of computing time spent on integrating the equations of motion is very small compared to the interaction computation usually involve extremely expensive double loops over all the particles in the system). Accuracy for large time steps is far more important, because the longer the time step that we can use, the fewer evaluations of the forces are needed per unit of simulation time.
\subsection{Centre of mass motion}
Centre of mass velocities are integrated through the following formula: \begin{equation} \mathbf{v}(t+\delta t/2)=\mathbf{v}(t-\delta t/2)+\frac{\delta t}{m} \mathbf{f}(t) \end{equation} 
 Coordinates are advanced through: \begin{equation} \mathbf{r}(t+\delta t)=\mathbf{r}(t)+\delta t\mathbf{v}(t+\delta t/2) \end{equation} 
The above equations are applied to each molecule~\cite[p~78-82]{allen}.

\paragraph{COM translation removal}
Being no external force acting on the system, the COM velocity should remain constant, say, zero. However, the integrator invariably develops a very slow change in the COM velocity. If such behaviour is not quenched, an appreciable COM motion develops eventually in long runs, and the temperature will be significantly misinterpreted\footnote{The same may happen due to overall rotational motion, but only when an isolated cluster is simulated. In periodic systems with filled boxes the overall rotation is coupled to other degrees of freedom and does not give any problems~\cite{gmx31}}. Hence in \brahms the COM net velocity is removed periodically.


The most popular integrator for solving Equation~\ref{eq:newton} is the venerable algorithm invented by Verlet:~\cite{verlet67} equivalent formulations are available, yielding identical trajectories but computing velocities with slightly improved accuracy. The algorithm employed in this work is the so-called {\em velocity Verlet} scheme:\cite{swope82a}
\begin{eqnarray} 
\mathbf{v}(t+\Delta t/2)&=&\mathbf{v}(t)+\Delta t\, \mathbf{f}(t)/2m \\ 
\mathbf{r}(t+\Delta t)&=&\mathbf{r}(t)+\Delta t\,\mathbf{v}(t+\Delta t/2)\\ 
\mathbf{v}(t+\Delta t)&=&\mathbf{v}(t+\Delta t/2)+\Delta t\, \mathbf{f}(t+\Delta t)/2m
\end{eqnarray} 
The above equations are applied to the mass centre of each site. Verlet-like schemes are employed in the great majority of MD programs due to their simplicity and exceptional stability; these algorithms are both symplectic and time-reversible%\footnote{It is impressive  that an almost 40 years-old numerical scheme is still state-of-the-art. It was gradually realized during the 19th century that the laws of motion embody extra, hidden, properties called {\em symplectic conservation}, in addition to more well-known conservation laws like conservation of energy and angular momentum. By the late 1960s the study of this structure gave rise to a new branch of geometry, symplectic geometry, which eventually, by the late 1980s, began to influence numerical analysts who developed special methods that obey the law of symplectic conservation. The Verlet method is the simplest method of this class, and still one of the best. At last, after being in popular use for several decades, the surprisingly good performance of the leapfrog method was explained. It's been possible to show that, although you can't stop the quantitative errors building up during the simulation, by obeying the law of symplectic conservation, you won't commit any qualitative errors and so the overall results will be far more reliable.}
.  


Equation~\ref{eq:euler} is more problematic than Newton's, since the treatment of rotational degrees of freedom is much less straightforward than that of linear motion; indeed, good integrators for Euler's equation are a relatively recent invention. 
One of the unique features of \brahms is the ability to handle efficiently rigid-body dynamics. Rigid bodies are not included in most simulation softwares because of the algorithmic complexity involved in propagating orientational degrees of freedom, as opposed to the relatively trivial problem of integrating the motion of point masses. Integrators that propagate rigid body orientational motion with an acceptable level of energy conservation for molecular dynamics are relatively new inventions.~\cite{meine05a}
Rigid body models possess symplectic structure  and time-reversal symmetry: standard numerical integration methods for rotational motion destroy both properties, introducing nonphysical dynamical behaviour and drift in the energy during long term simulations. For instance, the traditional and simple scheme by~\cite{finch81} (also reported in~\cite[p.~89]{allen}), based on a representation of rotations through quaternions\footnote{For a definition of quaternion and some notions on quaternion algebra, see Appendix~\ref{app:nlrb}.}, is neither symplectic nor reversible, and has indeed proved to yield a rather poor performance in terms of energy conservation, thus severely limiting the integration step size. 

In the last ten years, the mathematical community has developed more sophisticated but far superior rotational algorithms, that allow to  enlarge markedly the time-step, the key property being symplecticity\footnote{The development of symplectic rigid-body integrators is a significant methodology advance that generally regards non-symmetric rigid particles. The particular case of symmetric bodies can be readily treated with minor modifications, therefore it would be  possible to greatly enlarge the timestep in the simulation of Gay-Berne systems. However, symplectic integration appears to have gone unappreciated by the Gay-Berne community at large.}.
In particular, \cite*{dullw97a} have  proposed an integrator, that we shall call DLM,  based on the propagation of the full rotation matrix associated with each rigid body: in every time step, the orientational dynamics is integrated by a sequence of five planar rotations. The use of rotation matrices offers improved numerical stability, and since the method belongs to the leapfrog family of integrators, it means that  simple leapfrog integration techniques can be used for the entire set of dynamical equations appearing in the problem. 


The DLM integrator is symplectic and time-reversible, and conserves both the total linear and  angular momentum. %Furthermore, if the rigid-body is linear (e.g. a Gay-Berne site) then the five rotations reduce to one, thus lowering the computing time.
 The detailed algorithm is reported in Appendix~\ref{app:brahms}; it is interesting to note that the structure of the DLM scheme resembles the velocity Verlet algorithm, and indeed the two are identical in the way they treat linear motion.  Simulations on a system of 1000 rigid water molecules have shown that the step size allowed by the DLM scheme is about ten times larger than that possible with Fincham's scheme~\cite{fenne04a}. \cite{dullw97a} have also reported that the gap between their method and non-symplectic schemes increases with system size. The scheme has also been employed to simulate polymer helix formation using rigid-link methods~\cite{rapa02}.

%Other recent studies on rigid-body integration are the ones by~\cite{matub99a} (development of a reversible, though non-symplectic, integrator), \cite{mille02a} (reversible and symplectic scheme), and~\cite{kambe05a} (reversible and symplectic, but tested on linear molecules only). %It is interesting to note that all these algorithms rely on a quaternion representation, just as old~\cite{finch81}.

\paragraph{Gaussian themrostat Algorithm (velocity-Verlet form)} 
At first, the momenta of all molecules are advanced from $t$ to $t+\Delta t$:
\begin{equation} 
\label{eq:splitFirstStep1}
\mathbf{v}(t+\Delta t/2)=\mathbf{v}(t)+\Delta t\mathbf{f}(t)/2m \end{equation} 
\begin{equation} 
\label{eq:splitFirstStep2}
\mathbf{h}^b(t+\Delta t/2)=\mathbf{h}^b(t)+\Delta t \mathbf{T}^b(t)/2 \end{equation} 
where $\mathbf{h}^b$ and $\mathbf{T}^b=\mathbf{Q}(t)\mathbf{T}^S(t)$ are the body-frame angular momentum and torque, respectively.\\
Then mass centre positions and orientation matrices are moved a full time step according to the DLM scheme; subsequently, forces and torques at $t+\Delta t$ are computed.
Now the momenta at  $t+\Delta t$ are estimated:  
\begin{equation} 
\label{eq:splitFirstStep1}
\mathbf{v}(t+\Delta t)=\mathbf{v}(t+\Delta t/2)+\frac{\Delta t}{2} \mathbf{f}(t + \Delta t)/m \end{equation} 
\begin{equation} 
\label{eq:splitFirstStep2}
\mathbf{h}^b(t+\Delta t)=\mathbf{h}^b(t+\Delta t/2)+\frac{\Delta t}{2} \mathbf{T}^b(t+\Delta t) \end{equation} 
Then the Lagrange multiplier $\alpha$ is computed according to~\ref{eq:alpha}.


Jensen~\cite{jensen04} demonstrated that the CHARMM27 parameter set does not reproduce the experimental value for the area per lipid or the experimental order parameter profile when the bilayer is allowed to adjust its area freely, as in the N-P-T ensemble. In the N-P-T ensemble the area per lipid is severely underestimated and is steadily decreasing in time. Therefore, many membrane simulations using the CHARMM27 parameters are carried out in the N-P-A-T (A is the area of the bilayer in the xy plane) or N-Pz-$\gamma$-T ensembles~\cite{feller99}.
However, \cite{sonne05} have chosen to conduct all simulations in the N-P-T ensemble, using a new parameter set specifically developed.
\cite{benz05} simulated DOPC in the NPT ensemble, in a flexible-orthorombic cell, both with CHARMM27 and GROMACS. The CHARMM27 simulation resulted in an average d-spacing that was too high and an area/lipid that was (again) too low, suggestive of a bilayer that is less fluid than observed experimentally. The GROMACS simulation, on the other hand, yielded values that agreed reasonably well with experiment within experimental uncertainties.

\cite{ander80a} was the first to modify Newtonian equations of motion to sample ensembles other than the microcanonical: the volume is a dynamical variable while the generalised force acting on this variable is proportional to the difference between the internal and an externally fixed pressure. 


\subsection{Nose'-Hoover}
\cite{melchionna93} have extended the NPT scheme proposed by~\cite{hoove85a} to take into account a (size- and) shape-varying box. 
\subsubsection{Equations of motion}
Setting the external parameters $T_\mathrm{ext}$ and $P_\mathrm{ext}$ to the desired values of temperature and pressure, the isothermal-isobaric equations of motion can be written:  
\begin{eqnarray}
\dot{\mathbf{r}}_i&=&\mathbf{v}_i+\boldeta\,(\mathbf{r}_i-\mathbf{r}_\mathrm{COM})\\
\dot{\mathbf{v}}_i&=&\mathbf{f}_i/m_i-(\boldeta+\zeta\mathbf{1})\mathbf{v}_i\\
\dot{\mathbf{h}}^b_i&=&\mathbf{T}^b_i/I_i-\zeta\mathbf{h}^b_i\\% not sure
\dot{\zeta}&=&\left[T(t)/T_\mathrm{ext}-1\right]/\tau_T^2\\
\dot{\boldeta}&=&V\left[\mathbf{P}(t)-P_\mathrm{ext}\mathbf{1}\right]/(\mathrm{DOFs}\,k_B\,T_\mathrm{ext}\,\tau_P^2)\\
\dot{\mathbf{H}}&=&\boldeta\mathbf{H}
\end{eqnarray}
being:
\begin{itemize}
%\item $\eta$ tHe piston variable tHat acts as a strain-rate factor in order to balance the deviations of the instantaneous pressure from the preset value $P_\mathrm{ext}$.
\item $\boldeta$ the barostat, or rather a ``piston'' tensor variable that acts as a strain-rate factor in order to balance the deviations of the instantaneous pressure from the preset value~$P_\mathrm{ext}$ - the diagonal elements of~$\boldeta$ tend to adjust the internal pressure to the preset hydrostatic value~$P_\mathrm{ext}$, while the off-diagonal elements tend to eliminate nonzero fluctuations of the off-diagonal elements of the stress tensor~$\mathbf{P}$;
\item $\mathbf{r}_\mathrm{COM}=\sum_im_i\mathbf{r}_i/\sum_im_i$ the system centre of mass;
\item $\zeta$ the thermostat, or rather a ``frictional'' coefficient which evolves in time following the deviation of instantaneous temperature from the corresponding external parameter;
\item $\mathbf{1}$ is the identiy matrix\footnote{Explicitily, \begin{displaymath}
\mathbf{1}=\left(
\begin{array}{ccc}
 1 & 0 &0 \\
 0 & 1 & 0 \\
 0 & 0& 1
\end{array}
\right)
\end{displaymath} };
\item DOFs the total number of degrees of freedom of the system;
\item $T(t)$ the instantaneous temperature;
\item $\tau_T$ the time constant for relaxation of the temperature to the desired value - typically $\tau_T=1\div10\,$ps;
\item $\tau_P$ the time constant for relaxation of the pressure to the desired value - typically $\tau_P=10\div100\,$ps;
%\item $P(t)$ the sum of the kinetic and virial pressure terms;
\item $\mathbf{P}(t)$ the stress tensor;
\item $f$ the number of degrees of freedom in the system;
\item $\mathbf{H}=(\mathbf{a}, \mathbf{b}, \mathbf{c})$ the box matrix, where $\mathbf{a}, \mathbf{b}, \mathbf{c}$ are the basis vectors of the box.
\item $V$ the volume, given by\footnote{In general, the determinant of a $3\times3$ matrix is: $\det{[a_x b_x c_x; a_y b_y c_y; a_z b_z c_z]}= a_xb_yc_z-a_xc_yb_z-b_xa_yc_z+b_xc_ya_z+c_xa_yb_z-c_xb_ya_z$.} $V=\mathbf{c}\times\mathbf{b}\cdot\mathbf{a}=\det{\mathbf{H}}$
\end{itemize}
The values of $\tau_T$ and $\tau_P$ must be determined empirically: they have no particular physical meaning and are simply part of the computational technique. In principle, their values do not affect the final equilibrium results, but they do influence their accuracy and reliability, because if the fluctuations of kinetic energy and pressure are allowed to become too large it is difficult to think of the system existing in equilibrium at a particular temperature and stress~\cite{rapa}. 
The trace\footnote{The trace of a matrix is just the sum of the diagonal elements - more rigorously, the trace of an $n\times n$ matrix $A$ is defined as: $\mathrm{tr}(A)=\sum_{i=1}^na_{ii}$.} of $\boldeta=\mathbf{H}\mathbf{H}^{-1}$ has the property:
\begin{equation}
\mathrm{tr}(\boldeta)=\dot{V}/V 
\end{equation}
The pressure tensor is computed as: 
\begin{eqnarray}
\mathbf{P}(t)&=&\left[\sum_im_i\mathbf{v}_i(t)\otimes\mathbf{v}_i(t)\right]/V(t) + \mathbf{W}(t) \\
\mathbf{W}(t)&=&\left[\sum_i\sum_{j>i}\mathbf{r}_{ij}(t)\otimes\mathbf{f}_{ij}(t)\right]/V(t) 
\end{eqnarray}
with $\mathbf{W}(t)$ the stress tensor\footnote{Explicitely:
\begin{displaymath}
\mathbf{W}=\frac{1}{V}\sum_i\sum_{j>i}\left(
\begin{array}{ccc}
 r_{ij}^xf_{ij}^x & r_{ij}^xf_{ij}^y  &r_{ij}^xf_{ij}^z  \\
  r_{ij}^yf_{ij}^x&r_{ij}^yf_{ij}^y  &r_{ij}^yf_{ij}^z  \\
  r_{ij}^zf_{ij}^x& r_{ij}^zf_{ij}^y&r_{ij}^zf_{ij}^z 
\end{array}
\right)
\end{displaymath}}.
The instantaneous pressure is $P(t)=\mathrm{Tr}[\,\mathbf{P}(t)\,]\,/\,3$.

\subsubsection{Integration scheme}
The DLM scheme is extended as follows.
\paragraph{Part A}
\begin{eqnarray}
T(t)&:=&2\mathcal{K}[\mathbf{v}_i(t), \mathbf{h}_i^b(t)]\,/\,(k_B\,\mathrm{DOFs})\\
\mathbf{P}(t) &:=& \left[\sum_im_i\mathbf{v}_i(t)\otimes\mathbf{v}_i(t) + \sum_i\sum_{j>i}\mathbf{r}_{ij}(t)\otimes\mathbf{f}_{ij}(t)\right]\,/\,V(t)  \\
\mathbf{v}_i(t+\Delta t/2)&:=&\mathbf{v}_i(t)+\Delta t\,\left\{\,\mathbf{f}_i(t)/m\,-\,\left[\,\boldeta(t)+\zeta(t)\mathbf{1}\,\right]\cdot\mathbf{v}_i(t)\,\right\}\,/\,2\\
\mathbf{h}_i^b(t+\Delta t/2)&:=&\mathbf{h}_i^b(t)+\Delta t\,\left\{\,\mathbf{T}_i^b(t)\,-\,\left[\,\zeta(t)\mathbf{h}_i^b(t)\,\right]\right\}\,/\,2\\
\mathbf{Q}_i(t+\Delta t)&:=&\mathbf{Q}_i(t)\mathbf{R}_i^T[\mathbf{h}_i^b(t+\Delta t/2)]\\
\boldeta(t+\Delta t/2)&:=&\boldeta(t)+\Delta t\,V(t)\,\left[\,\mathbf{P}(t)-P_\mathrm{ext}\mathbf{1}\,\right]\,/\,\left[\,2Nk_BT(t)\tau_\mathrm{P}^2\,\right]\\
\mathbf{r}_i(t+\Delta t)&:=&\mathbf{r}_i(t)+\Delta t\,\left\{\,\mathbf{v}_i(t+\Delta t/2)\,+\,\boldeta\,(\,t+\Delta t/2)\,\cdot\,[\mathbf{r}_i(t+\Delta t)-\mathbf{r}_\mathrm{COM}]\right\}\\ 
\zeta(t+\Delta t/2)&:=&\zeta(t)+\,\Delta t\,[\,T(t)\,/\,T_\mathrm{ext}\,-\,1\,]\,/\,(2\tau_T^2)\\
\mathbf{H}(t+\Delta t)&:=&\mathbf{H}(t)\,\cdot\,\exp{[\Delta t\, \boldeta\,(t+\Delta t/2)]} \end{eqnarray}
Now forces and torques are computed at $t+\Delta t$; note that the volume is updated as $V(t+\Delta t)=\det{\mathbf{H}(t+\Delta t)}$.
\paragraph{Part B}
\begin{eqnarray}
T(t+\Delta t)&:=&2\mathcal{K}[\mathbf{v}_i(t+\Delta t), \mathbf{h}_i^b(t+\Delta t)]\,/\,(k_B\,\mathrm{DOFs})\\
\mathbf{P}(t+\Delta t) &:=& \left[\sum_im_i\mathbf{v}_i(t+\Delta t)\otimes\mathbf{v}_i(t+\Delta t) + \sum_i\sum_{j>i}\mathbf{r}_{ij}(t+\Delta t)\otimes\mathbf{f}_{ij}(t+\Delta t)\right]\,/\,V(t+\Delta t)  \\
\zeta(t+\Delta t)&:=&\zeta(t+\Delta t/2)+\,\Delta t\,[\,T(t+\Delta t)\,/\,T_\mathrm{ext}\,-\,1\,]\,/\,(2\tau_T^2)\\
\boldeta(t+\Delta t)&:=&\boldeta(t+\Delta t/2)+\Delta t\,V(t+\Delta t)\,\left[\,\mathbf{P}(t+\Delta t)-P_\mathrm{ext}\mathbf{1}\,\right]\,/\,\left[\,2Nk_B\,T(t+\Delta t)\,\tau_\mathrm{P}^2\,\right]\\
\mathbf{v}_i(t+\Delta t)&:=&\mathbf{v}_i(t+\Delta t/2)+\Delta t\,\left\{\,\mathbf{f}_i(t+\Delta t)/m\,-\,\left[\,\boldeta(t+\Delta t)+\zeta(t+\Delta t)\mathbf{1}\,\right]\cdot\mathbf{v}_i(t+\Delta t)\,\right\}\,/\,2\\
\mathbf{h}_i^b(t+\Delta t)&:=&\mathbf{h}_i^b(t+\Delta t/2)+\Delta t\,\left\{\,\mathbf{T}_i^b(t+\Delta t)\,-\,\left[\,\zeta(t+\Delta t)\,\mathbf{h}_i^b(t+\Delta t)\,\right]\right\}\,/\,2
\end{eqnarray}

\section{Radial Distribution Function}%A/T, p.54 - rapa p. 90 
The fluid state is characterised by the absence of any permanent structure. There are, nevertheless, well-defined structural correlations that can be measured experimentally to provide important details about the average molecular organisation. The radial (or pair) distribution function $g(r)$ - RDF for short - gives the probability of finding a pair of atoms a distance $r$ apart, relative to the probability expected for a completely random distribution at the same density. %; it is most simply thought of as the number of atoms a distance $r$ from a given atom compared with the number at the same distance in an ideal gas at the same density\footnote{Alternatively: the RDF is the ratio between the average number density at a distance $r$ from any given atom and the density at the same distance in an ideal gas at the same overall density.}. 
 The RDF can be evaluated through: \begin{equation} g(r)=\frac{V}{N^2}\left\langle\sum_i\sum_{j\ne i}\delta(r-r_{ij})\right\rangle=\frac{2V}{N^2}\left\langle\sum_{i}\sum_{j>i}\delta(r-r_{ij})\right\rangle \end{equation} 
where $V$ is the volume of the simulation region and $\delta$ is the Dirac function.
%The radial distribution function $g(r)$  describes the spherically averaged local organisation around any given atom. 

\section{Diffusion}%p122 rapa - p60 A/T 
Transport coefficients describe the material properties of a fluid within the framework of continuum fluid dynamics. Diffusion is the process whereby an initially nonuniform concentration profile (e.g. an ink drop in water) is smoothed in the absence of flow (no stirring). Diffusion is caused by the molecular motion of the particles in the fluid. The diffusion coefficient $D$, one of the most familiar transport coefficients, can be calculated from the Einstein expression (valid at long times): \begin{equation} 2tD = \frac{1}{3}\left\langle|\mathbf{r}_j(t)-\mathbf{r}_j(0)|^2\right\rangle \end{equation} where $\mathbf{r}_j(t)$ is the particle position. In practice, the average is computed for each of the $N$ particles in the simulation, the results added together and divided by $N$, to improve statistical accuracy. 
%: \begin{equation} D = \lim_{t \rightarrow \infty } \frac{1}{6Nt}\left\langle\sum_{j=1}^{N}\left[\mathbf{r}_j(t)-\mathbf{r}_j(0)\right]^2\right\rangle \end{equation} 
\subsection{Lipid lateral diffusion} 
In the literature, different units are
used. We use nm$^2$/$\mu$s as it makes for concise writing (no need
for exponentials) and intuitive meaning (nm and $\mu$s are about the
scales we investigate in the simulation). Conversions: $1\,\mathrm
{nm}^2/\mu \mathrm{s} = 1\,\mathrm{\mu m}^2/ \mathrm{s} = 10^{-8}\,\mathrm{cm}^2/\mathrm{s} = 10^{-12}\, \mathrm{m}^2/\mathrm{s} $.

\section{Dipole moment of a pair of equal and opposite point charges}
Given two equal yet opposite charges $q$, at a distance $l$ apart, the dipole moment is computed as:$$
\mu=ql
$$
The lipid headgroup dipole can be simply computed  in this way; the average is obtained by adding up and dividing the values  of all lipids.
 
\subsection{Head-group tilt angle}
Considering $\mathbf{d}$ the headgroup PN vector, the angle \dots

\subsection{Order parameters}
To describe the ordering of the hydrophobic lipid tails, the second order Legendre polynomial $P_2$ can be used:
\[P_2=\frac{1}{2}(3\cos^2\!\theta-1)\,,\]
where $\theta$ is the angle between the direction of the Gay-Berne tail site main axes $\mathbf{e}$ and the bilayer normal $\mathbf{n}$.
Since both $\mathbf{e}$ and $\mathbf{n}$ are unit vecotrs, we obtain $\cos^2\!\theta=(\mathbf{e}\cdot\mathbf{n})^2$.
Second-rank Gay-Berne order parameters $\overline{P}_{2}^{\,GB}$ have been calculated for the GB particles in the hydrocarbon chains.
 and compared to the experimental S$_{mol}$
segmental order parameter, which corresponds to the C-C bond order
parameter: $S_{CC}$. The Gay-Berne order parameter can be
defined as
$\overline{P}_{2}^{\,GB}=\left<P_{2}(\text{cos}\beta)\right>$, where
$\beta$ is the angle between the molecular axis and the director
\textbf{n} and $\left<\right>$ denotes a time average. 
Whithead et al. evaluated layer order parameter and analysed them versus temperature; it was clearly shown a phase transition from the solid (reminiscent of $L_{\beta'}$) to liquid-crystalline (reminiscent of $L_{\alpha}$) phase between $T^*=0.35$ and $T*=0.4$ (i.e. $290\,K<T<340\,K$), consistent with that identified previously from the diffusion coefficient data~\cite{white01a}.

\subsection{Elastic properties}

\subsubsection{Area compressibility modulus}
In the NPT ensemble (constant pressure and temperature), it is possible to estimate the area compressibility modulus~$K_A$ from the fluctuations in the lipid area.
Having computed the mean area per lipid $\langle A\rangle$ and the mean squared fluctuation  $\langle (A-\langle A\rangle)^2\rangle$ the compressibility modulus $K_A$ can be calculated as~\cite{feller99}
\begin{equation}
K_A=\frac{k_BT\langle A\rangle}{N_{leaflet}\langle (A-\langle A\rangle)^2\rangle}=\frac{k_BT\langle A\rangle}{N_{leaflet}\sigma^2(A)}
\end{equation}
with $N_{leaflet}$ the number of lipids per monolayer. 
This formula is widely used also when the pressure is controlled using the weak-coupling method~\cite{ber84}, which does not produce a perfect NPT ensemble. The effect of sampling from this spurious ensemble is negligible in this case, as:
\begin{itemize}
\item  the limitations set by small size and time scales are likely to be much more severe than this ensemble inaccuracy; 
\item   test simulations by~\cite{mar01} involving lipid bilayers using Parrinello-Rahman and Nose'-Hoover coupling schemes have not resulted in any noticeable effect on the observed fluctuations for correlation times larger than the coupling time, which is typically set to~$\sim0.5\,$ps. 
\end{itemize}
Considering area fluctuations, that occur on multi-nanosecond time scales, it is reasonable to assume that the weak-coupling method produces an isobaric ensemble.  
\begin{table}[]
\begin{center}
\vspace{8pt}
\begin{tabular}{lccc} % put @{} if you want to eliminate horizontal space between columns
&$\kappa$\,[bar$^{-1}$]&$\tau_P$\,[ps]& effective coupling strength\,[bar$^{-1}$ps$^{-1}$]\\
Bond et al. [2006] & $1\cdot10^{-5}$&40&$2.5\cdot10^{-7}$\\
\cite{bond06} & $1\cdot10^{-5}$&40&$2.5\cdot10^{-7}$\\
\cite{dicke05a} & $5\cdot10^{-6}$&1&$5\cdot10^{-6}$\\
\cite{falle04a} & $5\cdot10^{-6}$&1&$5\cdot10^{-6}$\\
\cite{shiqv05a} & $5\cdot10^{-6}$&1&$5\cdot10^{-6}$\\\hline
\end{tabular}
%\caption[Electron density comparison]{Total electron density profile: quantitative comparison. The head-head thickness $D_{HH}$ and selected electron densities from our simulations and experiments are compared.}
\label{tab:edp}
\end{center}
\end{table}

%\section{Monte-Carlo method} Along with MD, another important simulation tool is the {\em Monte-Carlo} (MC) method: here random steps are taken in order to achieve a rapid sampling of the most likely states of the system studied.  It should be stressed that only MD and not MC methods allow the theoretical possibility of obtaining time-dependent quantities from simulation, while both schemes can in principle be used for the same statistical-mechanical calculations.

\subsubsection{Volume compressibility modulus}
$K_V$ and can be related to the variance of the lipid volume $\sigma^2_{V}$:
\begin{equation}
K_V=\frac{k_BTV}{N_l\sigma^2_{V}} 
\end{equation}
with $V$ the average volume per lipid, $N_l$ the number of lipid molecules per layer, $T$ the temperature and $k_B$ the Boltzmann constant~\cite{venable06}.

\subsection{Lateral pressure distribution} Computer simulations of lipid bilayers indicate that bilayers have regions with negative lateral pressure trying to minimize the interfacial area, and regions with positive lateral pressure trying to expand the bilayer~\cite{goetz98a, lindahl00b, gulli04}. The mentioned MD bilayer studies, as well as less detailed models, predict lateral pressure variations in these regions of several hundred bars~\cite{sonne05}.

For inhomogeneous systems comprising an interface, hydrodynamics equilibrium requires that the component of the pressure tensor normal to the interface is constant throughout the system. The components tangential to the interface can vary in the interfacial region, but must be equal to the normal component in the bulk liquids~\cite[Sec.~17.1.2]{frenkel}.

For an inhomogeneous fluid there is no unambigous method to compute the pressure tensor components as a function of the interface normal. The popular Kirkwood-Buff convention involves dividing the system in $N$ equal slabs parallel to the interface, and weighing the virial over the slabs intersected by the line that connects every interacting pair of sites.



\section{Coda}
GBMOLDD~\cite{ilnyt01a}, DLPOLY~\cite{smith05c, todor04a, smith02b}, \cite{wilso97b, demig96a, wilson97} use non-symplectic integrators~\cite{fincham84, fincham93}.




\section{Controlling temperature and pressure}
Conventional MD simulations are carried out in the constant energy ensemble; on the other hand, it would be desirable to perform simulations in conditions closer to the real world, i.e. constant-Temperature and constant-Pressure.

The control can be achieved following two approaches:

\begin{itemize}
\item Constraint dynamics: T and P are kept rigorously fixed by constraints applied each step~\cite{evansHoover83}. 
\item Extended systems: additional degrees of freedom are introduced that couple dynamically to the rest of the system, thereby keeping T and/or P, {\em on average}, to the desired value~\cite{ander80a}.
\end{itemize}

 Pressure and temperature studies of biomembranes make possible the observation of new pressure-induced phases and changes in the dynamics of the lipid of which the membranes are composed~\cite{paci96}.
\subsubsection{Isotropic method}
In the isotropic control of pressure, the simulation volume retains its rectangular shape, so that changes consist of uniform contractions and expansions. % rapa, p 164
\subsubsection{Shape-changing method}
\cite{parrinello80} extended the case to a general simulation region in which the lengths and directions of the edges are allowed to vary independently, subject to uniform external pressure (an even more general case where the applied stress components are specified separately has also been developed~\cite{parrinello81}). 
The more flexible approach allows for the size and shape changes needed to accomodate lattice formation on freezing and for the study of structural phase transitions between different crystalline states~\cite[Section~6.2]{rapa}.

Conditions of constant pressure or stress allow the simulation of a solid-state phase transition with a change of cell size and shape. In a constant-stress simulation, the MD cell changes in size and shape in response to the imbalance between the internal and externally applied pressure~\cite{refson}.
Constraint methods~\cite{hoover82, evans83a}.

\cite{parrinello81, parrinello82} have developed MD methodology to apply to the system the most general condition of stress.


A metrix tensor $\mathbf{G}$ can be introduced:
\[\mathbf{G}=\mathbf{H}^T\mathbf{H}\]
Intersite distance can be computed:
\[\mathbf{r}_{ij}^2=\mathbf{s}^T_{ij}\mathbf{G}\mathbf{s}_{ij}\Rightarrow\mathbf{r}_{ij}=\sqrt{\mathbf{s}^T_{ij}\mathbf{G}\mathbf{s}_{ij}}\] 
where $\mathbf{s}_{ij}=\mathbf{s}_i-\mathbf{s}_j$  is the scaled vector specifying intersite distance in terms of lattice coordinates~\cite{hernandez01}.
The scaled coordinates span the unit cube and periodic images have coordinates $\mathbf{H}(\mathbf{s}+(n_x, n_y, n_z))$~\cite[Sec.~6.2]{rapa}.

To avoid overall cell rotation, the three sub-diagonal elements of $\mathbf{H}$ can be constrained to zero, that is, the $\mathbf{a}$ cell vector is constrained to lie along the $x$-axis and $\mathbf{b}$ is constrained to lie in the $xy-$plane~\cite{procacci97, refson}:
\begin{displaymath}
\mathbf{H}=\frac{1}{2}\left(
\begin{array}{ccc}
 a & b\,\cos\gamma & c\,\cos\beta \\
 0 & b\,\sin\gamma & c\,(\cos\alpha-\cos\beta\cos\gamma)/\sin\gamma \\
 0 & 0& (c/\sin\gamma)\sqrt{\sin^2\!\beta\,\sin^2\!\gamma-(\cos\alpha - \cos\beta\cos\gamma)^2}
\end{array}
\right)
\end{displaymath}
 Practically, at each time step the acceleration of those components, $\ddot{\mathbf{H}}_{ij}$, is set to zero (which is equivalent to adding a fictituous opposing force). This technique may also be used to allow uniaxial expansion only.


\subsection{Isotropic pressure tensor}
Controlling the pressure tensor so that $P_{xx}=P_{yy}=P_{zz}$, the area and volume can adjust and/or fluctuate separetely. This constraint, however, is equivalent to setting $\gamma=0$~\cite{feller95}.


\subsection{Changing box-shape}

Periodic boundaries and minimum image convention are most readily handled when the problem is expressed in terms of scaled coordinates, because the simulation region is then a {\em fixed unit cube}; use of physical variables introduces unnecessary complications when handling boundary crossings, because velocities and accelerations must be adjusted as well as coordinates.~\cite{rapa} %[Sec.~6.2, p.~158]
 It is then most convenient to introduce {\em scaled} (or {\em lattice}) coordinates $\mathbf{s}$ %, which give the position of atoms relative to the simulation cell
through the linear transformation:
\begin{equation}
\mathbf{r}=\mathbf{H}\mathbf{s}
\end{equation}
where $\mathbf{H}=(\mathbf{a}, \mathbf{b}, \mathbf{c})$ is a transformation matrix whose columns are the three vectors ($\mathbf{a}, \mathbf{b}, \mathbf{c}$) representing the edges of the simulation box. 
Conversly, a real-space vector $\mathbf{r}$ can be transformed into a box-space vector $\mathbf{s}$ via:
\begin{equation}
\mathbf{s}=\mathbf{H}^{-1}\mathbf{r}
\end{equation}
Since the cell vectors are linearly independent, matrix $\mathbf{H}$ can be inverted.


\subsection{Isotropic box}
The isotropic pressure $P_{iso}=\mathrm{Tr}[\,\mathbf{P}(t)\,]\,/\,3$ can be controlled by allowing {\em uniform} expansions and contractions of the system volume~\cite{ander80a}.
Obviously, in this case the ratio of surface area to system volume is implicitely fixed.

\section{Physical constants}

The physical constants employed by \brahms are reported in Table~\ref{tab:constants}. Some relations between units are collected in Table~\ref{tab:relations}.  % For the charge unit, note that $e=4.803\times10^{-10}$\,esu, and $e^2 = 331.8$\,(kcal/mol)\AA.

\begin{table}
\begin{center}
%\vspace{14pt}
\begin{tabular}{cll} % put @{} if you want to eliminate horizontal space between columns
Symbol & Name & Value  \\
$N_{AV}$ & Avogadro's number & $6.0221367\times10^{23}$\\
$k_{B}$ & Boltzmann's constant & $1.3806505\times 10^{-23}$\,J/K\\
$\varepsilon_0$ & Electric constant & $8.854 187 817\times 10^{-12}$\,F/m\\
$u$& Unified Atomic Mass Unit  &   $1.660 538 86\times10^{-27}$\,kg\\
$e$&elementary charge%\footnote{The elementary charge (symbol $e$ or sometimes $q$) is the electric charge carried by a single proton, or equivalently, the negative of the electric charge carried by a single electron. This is a fundamental physical constant and the unit of electric charge in the system of atomic units.}
&$1.602 176 53(14)\times10^{-19}$\,C
\end{tabular}
\caption[Physical constants]{Physical constants. Note: F/m = C$^2$N$^{-1}$m$^{-2}$, $e=4.803\times10^{-10}\,$esu, $e^2=331.8\,$kcal\,\AA/mole. % rapa, p 217
From~{\tt http://physics.nist.gov/cuu/Units/}.}
\label{tab:constants}  
\end{center}
\end{table}

\begin{table}
\begin{center}
%\vspace{14pt}
\begin{tabular}{cll} % put @{} if you want to eliminate horizontal space between columns
Symbol & Name & Relation to other units 
\\$V$ & Volt & $V = J / C$
\end{tabular}
\caption[Relations between units]{Relations between units.}
\label{tab:relations}  
\end{center}
\end{table}


\section{Configuration definition of the SSD water model} The atom coordinates of an SSD molecule located at $x=y=z=0$ and at orientation $\alpha=\beta=\gamma=0$ are reported in Table~\ref{tab:SSD coords}. 
\begin{table} \begin{center} \vspace{14pt} \begin{tabular}{lccc} % put @{} if you want to eliminate horizontal space between columns 
Atom & $x_1/$\AA &$x_2/$\AA&$x_3/$\AA\\ 
H$_1$& 0 & 0.75 & 0.53 \\ 
H$_2$& 0 & -0.75& 0.53\\ 
O &    0 & 0& -0.0654 
\end{tabular} 
\caption[SSD coordinates]{SSD coordinates in the principal (body-attached) frame of reference. Values from~\cite{bratk85a}.} 
\label{tab:SSD coords}  
\end{center} 
\end{table} 
The geometry of the SSD water is indeed equal to that of the TIP3P model~\cite{jorge81a}. Table~\ref{tab:PMI} reports the SSD principal moments of inertia corresponding to this geometry. 
\begin{table} \begin{center} \vspace{14pt} 
\begin{tabular}{lcc} % put @{} if you want to eliminate horizontal space between columns 
Moment of Inertia & Value in physical units & Reduced Value in \brahms 
\\ $I_y=m_Oz_0^2 + 2m_Hz_H^2$& 0.630\,amu\,\AA$^2$&0.003826\\ 
$I_z = 2m_Hy_H^2$&1.125\,amu\,\AA$^2$&0.006835\\ 
$I_x=m_Oz_0^2 + 2m_H(y_H^2+z_H^2)$&1.755\,amu\,\AA$^2$& 0.010661 
\end{tabular} 
\caption[SSD Principal Moments of Inertia in \brahms reduced units]{SSD Principal Moments of Inertia, expressed in reduced units. Of course, $I_x=I_y+I_z$.} 
\label{tab:PMI}   
\end{center}
\end{table}
For example, by setting $\sigma=3.4$\,\AA, $\epsilon/k_B=120$\,K and $m=39.95\times1.6747\times10^{-24}$\,g, then the MD time unit corresponds to $2.161\times10^{-12}$\,s.  \section{Dipolar potential} \begin{equation} \mathcal{V}_{ij}^{\mathrm{dp}}({\bf{r}}_{ij},{\Omega}_{i},{\Omega}_{j})= \frac{1}{4\pi\varepsilon\varepsilon_0}\left[ {\bmu_{i} \cdot \bmu_{j}}{r_{ij}^{3}} - \frac{3(\bmu_{i} \cdot {\bf{r}}_{ij})(\bmu_{j} \cdot {\bf{r}}_{ij})}{r_{ij}^{5}}\right]\;[=]\;\frac{Nm^2}{C^2}\frac{D^2}{m^3}\;[=]\;\frac{N}{C^2}\frac{C^2m^2}{m}\;\;[=]\; Nm \;[=]\; J \end{equation}




\subsection{RDF\index{RDF calculation}} The RDF computation has much in common with the interaction computation. The minimum image separations $r_{ij}$ of all the pairs of atoms are calculated and sorted into a histogram where each bin has a width $\Delta r$ and extends from $(n-1)\Delta r$ to $n\Delta r$, being $n$ the identifier of a given bin. If $h_n$ is the number of atom pairs ($i$, $j$) for which  $$(n-1)\Delta r \le r_{ij} < n\Delta r$$ then, assuming that $\Delta r$ is sufficiently small, we have the result: \begin{equation} g(r_n)=\frac{Vh_n}{2\pi N^2r_n^2\Delta r} \end{equation} where $r_n=(n-\frac{1}{2})\Delta r$. If the RDF measurements extend out to a maximum range $r_e$ the required number of histogram bins is $r_e/\Delta r$. The normalisation factors ensure that $g(r\rightarrow\infty)=1$. RDF functions are computed every $stepRdf$ steps; when $limitRdf$ computations have been performed and accumulated, results are averaged and outputted, and the calculation is re-initialised.   

\subsection{\index{Diffusion calculation}}
It is worth noting that in evaluating the diffusion coefficient it is important not to switch attention from one periodic image to another, which is why it is sometimes useful to have available a set of particle coordinates which have not been subjected to periodic boundary corrections during the simulation. For instance it is possible to keep track of the ``true'' atomic displacement $r_{j_x}'(t)$: \begin{equation}r_{j_x}'(t)=r_{j_x}(t)+nint((r_{j_x}'(t-\Delta t)-r_{j_x}(t))/L_x)\cdot L_x\end{equation}

The lateral diffusion coefficient can be computed 

\subsection{Quaternion-based integrators for rigid body dynamics}
An earlier implementation of \brahms employed quaternions to integrate rotations: this method showed poor stability, as it has been observed by other researchers~\cite{fenne04a, meine05a}.
In the following subsection the quaternion-based schemes are reported only for reference: they are no longer used in real simulations.
\subsubsection{Symmetric rigid bodies}
The leapfrog integration for the rotational motion of rigid linear sites needs at first the defining of a new entity, the ``gorque'' $\mathbf{g}$, as the equivalent force vector whose cross product with the orientation vector gives the torque\footnote{The torque $\mathbf{T}$ is therefore: $\mathbf{T}=\mathbf{e}\times\mathbf{g}$.}. Then we set $\mathbf{g}^\perp$ as:
\begin{equation} \mathbf{g}^\perp=\mathbf{g}-(\mathbf{g}\cdot\mathbf{e})\mathbf{e} \end{equation} 
the orientational velocity can be evaluated:
\begin{equation} \mathbf{u}(t+\delta t/2)=\mathbf{u}(t-\delta t/2)+\frac{\delta t}{I} \mathbf{g}^\perp(t) - 2[\mathbf{u}(t-\delta t/2)\cdot\mathbf{e}(t)]\mathbf{e}(t) \end{equation}
Orientations update:
\begin{equation}
\mathbf{e}(t+\delta t)=\mathbf{e}(t)+\delta t\mathbf{u}(t+\delta t/2)
\end{equation} 
The above equations are applied to each molecule~\cite[p.~90-91]{allen}.


\paragraph{NVT scheme} Constant translational-temperature dynamics is generated by the equations of motion: \begin{equation} \dot{\mathbf{r}}_i=\mathbf{p}_i/m_i\end{equation}\begin{equation} \dot{\mathbf{p}}_i=\mathbf{f}_i+\alpha\mathbf{p}_i \end{equation} \begin{equation} \label{lincem}\ddot{\mathbf{r}}_i=\frac{\mathbf{f}_i}{m_i}+\alpha\dot{\mathbf{r}}_i \end{equation} Imposing $\dot{K}=0$, or equivalently \begin{equation} \sum_im_i\dot{\mathbf{r}}_i\cdot\ddot{\mathbf{r}}_i=0\end{equation}  results in the following value for the Lagrange multiplier $\alpha$: \begin{equation} \alpha=-\frac{\sum_i\dot{\mathbf{r}}_i\cdot\mathbf{f}_i}{\sum_im_i|\dot{\mathbf{r}}_i|^2} \end{equation}\begin{equation} \alpha=-\frac{\sum_i\mathbf{p}_i\cdot\mathbf{f}_i}{\sum_i|\mathbf{p}_i|^2} \end{equation} The algorithm implemented consists at first in the computation of a projected (unconstrained) velocity $\mathbf{v}'_i$ through the Verlet scheme: \begin{equation} \label{com-firstStep}\mathbf{v}'_i(t)=\mathbf{v}_i(t-\delta t/2)+\frac{\delta t}{2m_i}\mathbf{f}_i(t)\end{equation} then the scaling factor $\beta\;(=[1+\alpha\,\delta t/2m]^{-1})$ is computed as follows: \begin{equation} \label{com-secondStep}\beta = \sqrt{\frac{3(N-1)k_BT_{trans}}{\sum_{i=1}^Nm_i|\mathbf{v}'_i(t)|^2}}\end{equation} where $N$ is the {\emph total} number of sites. Eventually, the integration step is: \begin{equation} \label{com-thirdStep}\mathbf{v}_i(t+\delta t/2)=(2\beta - 1)\cdot\mathbf{v}_i(t-\delta t/2)+\beta\frac{\delta t}{m_i}\mathbf{f}_i(t)\end{equation} \paragraph{Rotational motion - Linear rigid sites}  Hoover's method can be extended further in order to take into account linear molecules, so to constraint the rotational kinetic energy as well. The constrained rotational equation of motion is: \begin{equation} \label{rotcem} \dot{\mathbf{\omega}}_i=\frac{\mathbf{\tau}_i}{I_i}+\alpha{\mathbf{\omega}}_i \end{equation} Combining Eq.~\ref{rotcem} with Eq.~\ref{lincem} and imposing $\dot{K}=0$, or equivalently \begin{equation} \sum_i ( m_i\dot{\mathbf{r}}_i\cdot\ddot{\mathbf{r}}_i + I_i{\mathbf{\omega}}_i\cdot\dot{\mathbf{\omega}}_i)=0\end{equation}  results in the following value for the Lagrange multiplier $\alpha$: \begin{equation} \alpha=-\frac{\sum_i(\dot{\mathbf{r}}_i\cdot\mathbf{f}_i +{\mathbf{\omega}}_i\cdot\mathbf{\tau}_i)}{\sum_i( m_i|\dot{\mathbf{r}}_i|^2 + I_i|\mathbf{\omega}_i|^2)} \end{equation} The constrained equation for integrating the rotational motion is:  \begin{equation}\dot{\mathbf{u}}_i=\mathbf{g}^\perp_i/I_i-\lambda \mathbf{e}_i+\alpha{\mathbf{u}}_i \end{equation} where $\mathbf{g}^\perp$ is the component of the equivalent force perpendicular to the orientation vector: \begin{equation}\mathbf{g}_i^\perp=\mathbf{g}_i-(\mathbf{g}_i\cdot\mathbf{e}_i)\mathbf{e}_i\end{equation} So in this case there are two constrains, $\alpha$ and $\lambda=\mathbf{u}^2$.  algorithm implemented consists at first in the computation of a projected (unconstrained) velocity $\mathbf{u}'_i$ through: \begin{equation} \label{rot-firstStep}\mathbf{u}'_i(t)=\mathbf{u}_i(t-\delta t/2)-[\mathbf{u}(t-\delta t/2)\cdot\mathbf{e}(t)]\mathbf{e}(t)+\frac{\delta t}{2I_i}\mathbf{g}^\perp_i(t)\end{equation} e scaling factor $\beta\;(=[1+\alpha\,\delta t/2]^{-1})$ is computed as follows: \begin{equation} \label{rot-secondStep}\beta = \sqrt{\frac{2NT_{rot}}{\sum_{i=1}^NI_i|\mathbf{u}'_i(t)|^2}}\end{equation} where $N$ is the {\emph total} number of sites. Eventually, the integration step is: \begin{equation} \label{rot-thirdStep}\mathbf{u}_i(t+\delta t/2)=(2\beta - 1)\cdot\mathbf{u}_i(t-\delta t/2)+\beta\frac{\delta t}{I_i}\mathbf{g}^\perp_i(t)-2\beta[\mathbf{u}(t-\delta t/2)\cdot\mathbf{e}(t)]\mathbf{e}(t)\end{equation}

\paragraph{NPT scheme} The NPH constant-pressure method by~\cite{ander80a} can be combined with the NVT scheme previously considered: the resulting NPT algorithm has been implemented in the code\footnote{The main reference employed is the paper by~\cite{brown94a}, where the equations of motion are presented in a readily implementable form.}.

\subparagraph{Volume V}
\begin{equation}
V(t+\delta t) = 2V(t) - V(t-\delta t) + \frac{\delta t^2}{Q}\left[P(t)-P_E\right]
\end{equation}
\begin{equation}
\dot{V}(t) = \frac{V(t+\delta t)-V( t-\delta t)}{2\delta t}
\end{equation}
\begin{equation}
\ddot{V}(t) = \frac{P-P_E}{Q}
\end{equation}
$P_E$ is the external pressure, $V$ is the volume of the system and $Q$ is a constant\footnote{The ``piston mass'' $Q$ is an adjustable parameter in Andersen's method; a low ``mass'' will result in rapid box size oscillations, which are not damped very efficiently by the random motions of the molecules. A large ``mass'' will give rise to slow exploration of volume-space, and an infinite mass restores normal MD~\cite[p.~234]{allen}.}.
$Q$ largely influences the rate of relaxation of the pressure and also the magnitude of the fluctuations of the pressure. As~$Q\rightarrow\infty$, the NPH ensemble tends towards the NVE ensemble. Static properties are relatively insensitive to the value of $Q$ within fairly broad limits; small values of $Q$ promotes faster equilibration to $P_E$, but can lead to nonphysically rapid fluctuations in the volume or to (explosion) catastrophic irreversible expansion of the system\footnote{In their simulation runs on a LJ system, \cite{brown94a} set $Q^{*}(=Q\sigma^4/m)=0.0027$.}. 
\subparagraph{Translational (centre of mass) motion}
Computation of a projected (unconstrained) velocity $\mathbf{v}'_i$ through the Verlet scheme:
\begin{equation} \label{npt-firstStep}\mathbf{v}'_i\left(t\right)=\mathbf{v}_i\left(t-\delta t/2\right)+\frac{\delta t}{2}\left[\frac{\mathbf{f}_i\left(t\right)}{m_i}+\frac{\mathbf{r}_i\left(t\right)}{3V\left(t\right)}\left(\ddot{V}\left(t\right)-\frac{2\dot{V}^2(t)}{3V(t)}\right)\right]\end{equation}
then the scaling factor $\beta$ is computed as follows:
\begin{equation} \label{npt-secondStep}\beta = \sqrt{\frac{3\left(N-1\right)k_BT_{trans}}{\sum_{i=1}^Nm_i|\mathbf{v}'_i\left(t\right)|^2}}\end{equation}
where $N$ is the {\emph total} number of sites. The leapfrog step for updating the velocities becomes:
\begin{equation} \label{npt-thirdStep}\mathbf{v}_i\left(t+\frac{\delta t}{2}\right)=\left(2\beta - 1\right)\mathbf{v}_i\left(t-\frac{\delta t}{2}\right)+\left(1-\beta\right)\mathbf{r}_i\left(t\right)\frac{2\dot{V}\left(t\right)}{3V\left(t\right)}+\beta\delta t\left[\frac{\mathbf{f}_i\left(t\right)}{m_i}+\frac{\mathbf{r}_i\left(t\right)}{3V\left(t\right)}\left(\ddot{V}\left(t\right)-\frac{2\dot{V}^2(t)}{3V(t)}\right)\right]\end{equation}



\subsubsection{General (non-symmetric) rigid bodies}
As it has long been common practise in MD, the orientations of rigid non-linear molecules can be represented by {\em quaternions} $\mathbf{q}$, ordered number quartets:~\cite{finch81} developed an algorithm that propagates the rotational dynamics via the update of the quaternion associated with each site.   
The first step is to bring all the angular momenta $\mathbf{h}$ up to date:
\begin{equation}
\mathbf{h}^S(t)=\mathbf{h}^S(t-\delta t/2)+\frac{\delta t}{2}\mathbf{T}^S(t)
\end{equation}
Then a guess at $\mathbf{q}(t+\delta t/2)$ is made:
\begin{equation}
\mathbf{q}(t+\delta t/2)=\mathbf{q}(t)+\frac{\delta t}{2}\dot{\mathbf{q}}(t)
\end{equation}
The main algorithm equations are~\cite[p.~89]{allen}:
\begin{equation}
\mathbf{h}^S(t+\delta t)=\mathbf{h}^S(t-\delta t/2)+\delta t\mathbf{T}^S(t)
\end{equation}
and:
\begin{equation}
\mathbf{q}(t+\delta t)=\mathbf{q}(t)+\delta t\dot{\mathbf{q}}(t+\delta t/2)
\end{equation}
At every step, all quaternions have to be adjusted in order to preserve the constraint $q_1^2+q_2^2+q_3^2+q_4^2=1$: this operation breaks time-symmetry, hence the algorithm is not reversible. Note that Fincham's algorithm is not even symplectic: indeed it has shown poor performances in terms of long-time energy conservation.

\paragraph{Twin-range method}
When simulating species with a significant electrostatic contribution, it may be desirable to use a twin-range method, in which two cutoffs are specified: all interactions below the lower cutoff are calculated as normal at each step, whereas the interactions due to atoms between the lower and upper cutoffs are evaluated only when the neighbour list is updated and are kept constant between these updates~\cite[Sec.~6.7.1]{leach}.
\paragraph{Group-based cutoff}
Using simple atom-based cutoffs the interaction energy may fluctuate violently near the cutoff distance. This problem can be avoided by calculating the interactions on a group-group basis, being the  ``group'' a collection atoms of zero total charge (e.g. all three atoms in a water molecule)~\cite[Sec.~6.7.2]{leach}. 

\paragraph{Long-range effects}
In all-atom simulations, the majority of the MD community  is currently oriented towards including long-range electrostatics via Ewald methods. Nonetheless, some researchers are openly {\it against} the use of Ewald methods~\cite{scott02a, beckd05a}.
For example, considering atomic-level membrane studies, some authors have tested both cutoff and Ewald sums and concluded that cutoff simulations do not suffer to any great extent from the omission of long-range electrostatics~\cite{marri01c,rog03,wohle04a}.
On the other hand, \cite{patra03a, patra04a} have shown that long-range truncation leads to artifacts such as enhanced order of acyl chains together with decreased diffusion and areas per lipid (with respect to experimental data). \cite{polya05a} also found that Ewald performs better than cutoff, but interestingly, their results in some cases showed opposite trends compared to Patra and coworkers.

As far as we are concerned, the decision has been made to avoid the expensive computation of long-range forces in favour of a faster cutoff scheme, this choice being also consistent with the overall coarse-grained spirit of the project.




\section{Methodological issues in lipid bilayer simulation}

In this section we consider  specific methodology aspects important to simulate bilayer systems. 

\subsection{Equilibration}\index{Equilibration of membrane - MD}

Recent MD studies of hydrated biomembranes have shown that some important bilayer properties (e.g. the lipid area) fluctuate on the time scale of several nanoseconds~\cite{lindahl00, mar01, anezo03a, hogberg06}. 
Therefore, to obtain reliable equilibrium properties, membrane systems must be simulated for realtively long times, at least a few tens of nanoseconds\footnote{In a time where 1-ns (and less) NPT trajectories were published without too much worries, \cite{feller99} profhetically warmed that ``the inherently slow process of area fluctuation requires simulations close to 10\,ns, rather than~1\,ns, before convergence and/or the generation of sufficient equilibrium sampling is obtained''.}. 

\paragraph{\index{Surface tension}}
The surface tension $\gamma$ along a bilayer is usually calculated from the difference between the normal ($P_N=P_{zz}$) and lateral ($P_L=(P_{xx}+P_{yy})/2$) components of the pressure tensor:
\begin{equation}
\gamma=\int[P_N(z)-P_L(z)]\,\de z
\end{equation}
If sufficient water is placed outside the bilayer, the integrand vanishes at the endpoints, since bulk water is tension-free: we can therefore integrate over the entire simulation volume without making any arbitrary definition of the interface. %~\cite{gulli04}.
The surface tension per leaflet of a bilayer 
is related to the pressure tensor by~\cite{feller95}:
\begin{equation}\label{eq:surfTens} 
\boxed{\gamma= \langle L_z\times(P_{zz}-P_{xx}/2-P_{yy}/2)\rangle\,/\,2}
\end{equation}
where $L_z$ and $P_{zz}$ denote the length of the unit cell and component of the pressure tensor normal to the surface, respectively, and $P_{xx}$ and $P_{yy}$ are tangential components of the pressure tensor. Positive values of the surface tension result from tangential pressures that are less than the pressure along the normal~\cite{feller99}. Note that eq.~\ref{eq:surfTens} gives the surface tension {\em per leaflet}; typical units are dyn/cm/leaflet~\cite{pastor05}. 


\subsection{Ideal conditions for membrane simulation}
Simulations have also been performed in the NP$_n$AT ensemble (constant mole number, pressure normal to the bilayer plane, area/lipid, temperature) in which the area/lipid is restrained about the experimentally determined value~\cite{bempo04a, bempo04b}. However, to perform such simulations, one must have prior experimental knowledge of the area/lipid of the system being studied. Complicating matters, this parameter has been shown to vary greatly with lipid type and hydration level~\cite{nagle00a}, so each system simulated at NP$_n$AT requires an area/lipid value specific for that system, which may not be available. Consequently, NP$_n$AT simulations will be of marginal value if simulations are to be used eventually in place of experimental analyses of new, unstudied bilayer systems, such as those containing mixtures of lipids and/or membrane proteins. Simulations performed in the NPT (constant mole number, pressure, and temperature) ensemble have the potential for accomplishing this objective, because, in principle, they do not require prior experimental information: this makes it possible to verify the area per lipid as a result of simulations of other lipids or lipid-protein systems~\cite{tiele97a}. Historically, a pressure of 1 \,bar (both lateral and perpendicular components) has been used: this is actually the same as zero bar within the numerical accuracy of pressure calculation. Given perfect force fields and constant pressure and temperature algorithms, an NPT simulation should be adequate for reproducing accurate experimental results~\cite{benz05}, so that the size and the shape of the simulation box are free to adjust, allowing the membrane area and thickness to fluctuate. This offers the possibility to compare and validate simulations by examining their ability to reproduce important structural quantities such as the projected area per lipid. 

Some workers control the pressure by separate coupling to 1\,atm in the three space coordinates (normal and lateral directions), corresponding to a stress-free bilayer; the simulation region is constrained to orthorombic symmetry~\cite{marri01c, anezo03a, wohle04a, benz05}.

In some MD simulations, DPPC bilayers of 36-40 lipids/leaflet near the experimental surface area per lipid have shown surface tensions near 20\,dyn/cm. This result is controversial because the experimental surface tension of macroscopic lipid bilayers is close to or equal to zero~\cite{jahnig96}. 
\cite{feller96} propose that $\gamma\ne0$ for simulation-sized bilayers because long-wavelength undulations are not present. Hence a nonzero surface tension must be applied in NP$\gamma$T simulations in what may be regarded as a small system correction.
Others have argued that the nonzero surface tension observed in simulations is an artifact for the force field, and hence the appropriate boundary condition for lipid bilayer simulations is constant isotropic pressure~\cite{jahnig96, tieleman96, berge97a, lindahl00b, marrink01, benz05, sonne05}; from this point of view, an ideal force field would yield $\gamma=0$\, dyn/cm for a non-stressed lipid bilayer at the correct surface area. 

\cite{feller95}

In the simulations by~\cite{anezo03a}, the pressure was controlled either anisotropically or semi-isotropically. In the former case, the three unit-cell dimensions fluctuate independently from each other, and the total pressure P remains constant. This corresponds to an NP$_x$P$_y$P$_z$T ensemble, which is not rigorously defined and stable only when at least two of the pressure components are equal. The semi-isotropic case corresponds to NP$_n$P$_l$T. In both cases, the pressure components are kept at 1\,bar on average. The only difference in the simulations is that in the anisotropic case the simulation box fluctuates independently in x and y whereas in the semi-isotropic case the interface maintains a square. When the lateral pressure and normal pressure are equal, the average surface tension is zero; at constant box length, specifying zero surface tension and a normal pressure of 1 bar is equivalent to specifying a lateral and normal pressure of 1 bar.  \cite{anezo03a} ignore the effect of the fluctuating box length and assume that specifying a lateral and normal pressure of 1 bar is the same as specifying zero surface tension and a normal pressure of 1 bar.

A comparison between the original and shifted-force form of the Coulomb potential is presented in Fig.~\ref{fig:sf}.
\begin{figure}
\centering
\mbox{\subfigure[]{\includegraphics[width=75mm]{methodology/shiftPot.eps}}\quad
      \subfigure[]{\includegraphics[width=75mm]{methodology/shiftForceInset.eps}}}
\caption[Coulomb interaction: normal and shifted-force]{Coulomb interaction. Comparison between the natural and the shifted-force forms. (a) Potential. (b) Force.
}
\label{fig:sf}
\end{figure}

%If appropriately constructed shift functions are used for the electrostatic forces, no charge groups are needed, but in principle, the cutoff scheme should be combined with a lattice sum for long-range electrostatics~\cite{gmx31}. \cite{brooks85} found that the shifted force scheme (Equation~\ref{eq:sf}) resulted in a faithful representation of short-ranged structures.

\subsubsection{Surface tension}\index{lateral pressure profile!surface tension}
The surface tension per leaflet of a bilayer 
is simply related to the pressure tensor by:~\cite{feller95}
\begin{equation}\label{eq:surfTens} 
\gamma= \langle L_z\times(P_{N}-P_{L})\rangle\,/\,2
\end{equation}
where $L_z$ and $P_{N}$ denote the length of the unit cell and component of the pressure tensor normal to the surface, respectively, and $P_{L}$ is the lateral components of the pressure tensor. 
The surface tension per monolayer $\gamma$ can be also computed from the lateral pressure profile as:
\begin{equation}
\gamma=\frac{1}{2}\int_0^{L_z}[P_{N}(z)-P_{L}(z)]\:\de z=-\frac{1}{2}\int_0^{L_z}\pi(z)\:\de z
\end{equation}
The surface tension measures the intrinsic tendency of the bilayer lateral dimensions to change relative to the longitudinal one, while the total pressure measures the tendency of the three dimensions to change simultaneously. A positive surface tension indicates that the system tends to a lateral shrinkage and a longitudinal expansion. A zero surface tension indicates that the system is at equilibrium with respect to its surface area and stacking distance for a particular pressure.~\cite{devri05a}



\subsection{Constraint method - Gaussian thermostat}
\label{ss:hoov} 
The Gaussian thermostat method~\cite{hoover82,evans83} enforces constant temperature by introducing %nonholonomic 
constraints into the equations of motion to fix the kinetic energy; in effect this serves as a mathematical thermostat.~\cite{rapa} The equilibrium properties of this isothermal system can be shown to be those of the canonical ensemble. %The constraint-temperature method has been implemented in the code; the main references employed are the papers by~\cite{brown94a} (for the centre of mass motion) and by~\cite{finch86a} (extension to take into account linear molecules and hence rotational motion). 
The total kinetic energy $\mathcal{K}$ is: \begin{equation} \mathcal{K}=\frac{1}{2}\sum_im_i|\mathbf{v}_i|^2+\frac{1}{2}\sum_i\sum_xI_{ix}(\omega^b_{ix})^2\end{equation}
 being, for each particle $i$, $m_i$ the mass, $\mathbf{v}_i$ the linear velocity, $I_{ix}$ the $x-$component of the principal moments of inertia,  $\omega^b_x$ the $x-$component of the body-frame angular velocity; $\sum_i$ denotes a sum over the system's particles and $\sum_x$ denotes a sum over the vector components. 

by setting $\dot{\mathcal{K}}=0$ we obtain the Lagrange multiplier $\alpha$:
\begin{equation}\label{eq:alpha}
\alpha = - \frac{\sum_i\mathbf{v}_i\cdot\mathbf{f}_i+\sum_i\bomega^b_i\cdot\mathbf{T}^b_i}{\sum_im_i\mathbf{v}^2_i+\sum_i\sum_xI_{ix}(\omega^b_{ix})^2}
\end{equation}
being $\bomega^b$ and $\mathbf{T}^b$ the body-frame angular velocity and torque; $\alpha$ has the dimension 1/time.  
Forces and torques are modified according to:
\begin{equation*}
\mathbf{f}_i := \mathbf{f}_i +\alpha\, m_i\mathbf{v}_i
\end{equation*}
\begin{equation*}
\mathbf{T}_i := \mathbf{T}_i +\alpha\,\mathbf{h}_i
\end{equation*}
The constraint method thus reduces to simple scaling of the forces as well as the velocities. %; for other integrators this may not be the case.
It must be noticed that the constraint method is only designed to conserve $\mathcal{T}$, not to fix this quantity at the desired value $T$ (which do not in fact appear in the algorithm). This means that $\mathcal{T}$ will tend to drift and it could be necessary to make small corrections every $n$ time-steps. %AllentTildesley 239: addirittura *every* time step! 
This is done through simple velocity scaling.   


%\subsection{Nose'-Hoover thermostat} The popular Nose'-Hoover algorithm is described in the next section, where it will be considered coupled to a barostat for pressure control. Temperature and pressure are not kept fixed by this scheme; instead, using a negative feedback, the fluctuations are limited and the mean values are equal to the desired values\footnote{This actually makes more physical sense, as in reality we do expect some fluctuations.}.


\paragraph{Symmetric rigid-bodies} The DLM method can be modified to specifically treat the symmetric rigid body, with (say) $I_1=I_2$ in principal coordinates. \begin{eqnarray*} \mathbf{R}_1:=\mathbf{R}_z\left(\;\Delta t\,\left(\frac{1}{I_3}-\frac{1}{I_2}\right)\;\right)&\;\;\rightarrow\;\; \mathbf{h}^b=&\mathbf{R}_1\mathbf{h}^b\\ \mathbf{R}_2:=\mathbf{R}_{\mathbf{h}^b}\left(\;\Delta t\,\frac{|\mathbf{h}^b|}{I_2}\;\right)&\;\;\rightarrow\;\; \mathbf{Q}^T=&\mathbf{Q}^T\mathbf{R}_2^T; \end{eqnarray*} 
These formulae allow the treatment of, for instance, point-dipoles and Gay-Berne solid ellipsoids. 
